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Abstract 

We present a rational approximation for rapid and accurate com¬ 
putation of the Voigt function, obtained by residue calculus. The 
computational test reveals that with only 16 summation terms this 
approximation provides average accuracy 10“^^ over a wide domain 
of practical interest 0 < x < 40,000 and 10“^ < y < 10^ for ap¬ 
plications using the HITRAN molecular spectroscopic database. The 
proposed rational approximation takes less than half the computation 
time of that required by Weideman’s rational approximation. Algo¬ 
rithmic stability is achieved due to absence of the poles at y ^ 0 and 
—OO < X < oo. 

Keywords: Voigt function, Faddeeva function, complex probability 
function, complex error function, rational approximation, spectral line 
broadening 


1 Introduction 

The Voigt function is widely used and finds broad applications in many scien¬ 
tific disciplines [H El El lU [5] . It is commonly applied in Applied Mathematics, 
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Physics, Chemistry and Astronomy as it describes the line prohle behavior 
that occnrs dne to simnltaneons Lorentz and Doppler broadening effects; 
the Lorentz broadening is observed as a resnlt of the Heisenberg nncertainty 
principle and chaotic mnltiple collisions of the particles while the Doppler 
broadening appears dne to velocity distribntion of the particles. 

The Voigt fnnction can describe the spectral properties in the photon 
emission or absorption of atmospheric gases 0 El El El fin] and celestial 
bodies m- It is also widely nsed in crystallography na and can be ntilized 
in many other spectroscopic applications, for example, to characterize the 
photo-lnminescent properties of nanomaterials [13] or to determine the hyper 
strnctnre of an isotope [T3| and so on. 

Mathematically, the Voigt fnnction is a convolntion integral of the Canchy 
and Ganssian distribntions [H El El HI E] 


K {x,y) 



- T^dt 

1/2 + (a; - t) 


( 1 ) 


and represents the real part of the complex probability fnnction EE] 



—oo 


where z = x + iy vs a. complex argument. The complex probability function 
can be expressed explicitly as a superposition of the real and imaginary parts 
W {x, y) = K (x, y) + iL {x, y), where its imaginary part is given by [21 E] 


L{x,y) 


1 7 

J 1/2 + (x — ty 

— OO 


(3) 


Another closely related function is the complex error function, also known 
as the Faddeeva function 0 da ESI Enm uni 


w{z) = e ^ [1 — erf {—iz)] 


= e " I 1 + f e^^dt 


(4) 
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There is a relation between complex probability function (§ and complex 
error function Q 

W {z) = w {z), Im [ 2 ;] ^ 0 

or 

w{x,y) = K {x,y)+iL{x,y), y > 0- (5) 

'-V-^ 

W{x,y) 

In order to describe spectral characteristics of a system with high resolu¬ 
tion, intense computation is required. For example, in a line-by-line radiative 
transfer simulation to resolve some problems associated with inhomogeneity, 
the Earth’s or other planetary atmosphere can be divided up to 1000 layers 
PE]. Taking into account that computation requires a nested loop procedure 
in order to adjust properly for the htting parameters for each atmospheric 
layer that may contain many different molecular species, the total number 
of the computed points may exceed hundreds of millions. Since in a radia¬ 
tive transfer model the computation of spectral broadening prohles requires 
considerable amount of time, a rapid approximation of the Voigt function is 
very desirable EE!. Consequently, the rapid and accurate computation of 
the Voigt/complex error function still remains topical (see for example an 
optimized algorithm in the recent work [2U]h 

In this work we present a new rational approximation of the Voigt func¬ 
tion for efficient computation. Due to absence of the poles at y ^ 0 and 
—oo<x<oo this rational approximation enables stability in algorithmic 
implementation. 

2 Derivation of the rational approximation 

The complex error function (|^ can also be expressed in an alternative form 
as [211122] 

00 

w {x, y) = J exp (—f^/4) exp {—yt) exp (ixt) dt, 

0 

where its real and imaginary parts are 

00 

Re [w {x, y)] = J exp (—f^/4) exp {—yt) cos (xt) dt 
0 
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and 


OO 


Im [w {x, y)] = J exp (—4) exp {—yt) sin (xt) dt, 

0 

respectively. By changing sign of the variable x to negative in the last two 
eqnations above, we can see the symmetric properties of the complex error 
fnnction 

f Re [w {x, y)] = Re [w {-x, y)] 

1 Im [w {—X, y)] = —Im [w {x, y)]. 

Conseqnently, it follows that 


Re [w (x, y)] = [w (x, y) + w (-x, y)] /2 (6) 

and 

Im [w (x, y)] = [w (x, y) - w (-x, y)] /2. 

It is worth noting that with these identities and eqnation Q we can also 
obtain two interesting relations for the real and imaginary parts of the error 
fnnction of complex argnment as follows 


{ Re [erf (x + iy)] 
Im [erf (x + iy)] 


erf (x + iy) + erf (x — iy) 
2 

erf (x + iy) — erf (x — iy) 
2i 


Since P 
where 6 {x — t) 


lim •- = nd (x — t) e * 

y^O 2/2 + (a; _ 

is the Dirac’s delta fnnction, we obtain 




y 


TT y y'^ + [x — ty 


jdt 


ye 




= lim — / 

j/^o n J 2/2 + (x - ty 


:dt 


y=0 


= — / ttS {x — t) e ^ dt = e 
vr ./ 


Conseqnently, we can write Re [w (x, y = 0)] = K {x,y = 0) = exp (—x^) and 
from the identity (|^ it immediately follows that 

exp (-x^) = [K {x,y = 0) + K (-x, y = 0)] /2. (7) 
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In our recent publication we have shown that a sampling methodology 
based on incomplete expansion of the sine function leads to a new series 
approximation of the complex error function [23] 


2 M -1 

W {z) =W {z) K, ^ 

m=l 


+ {z + k / 2 ) Bm 
Cl-{z + i<;/2f ’ 


Im [ 2 ;] ^ 0 . 


( 8 ) 


where the coefficients are 


^/n{2m-l) ^ /vr ( 2 m - 1 ) (nh + < 2 / 2 ) 


Am. — 


22 ^h 


n=—N 
N 


sm 


2^h 


Bm - 2 M -1 


1 , J 


( 7r(2m-l) (nh + .?/ 2 ) 


n=-N 


2^h 


Cm = 


71 {2m — 1 ) 


2M+ih 

with <2 = 2.75, h = 0.25, M = 5 and N = 23. As we can see, the integer 
on upper limit of the summation in this approximation is equal to 2^~^. 
However, this restriction can be omitted and application of the series ap¬ 
proximation above can be generalized for any integer. 

Consider the following limit for the sine function 


oM —1 oM —1 

1 /2m -1 

mAA 2 ^ ^ V 2 ^ 7 “ M “co ^ 

m=l ^ ^ m=l 


5inc(t)= lim ^ E 


COS 


m — 1/2 


\ 2 ^ ) M^oo 2^-1 V 2^ 

where we imply that the sine function is defined as 

{sine (t 7 ^ 0 ) = sin {t) /t, sine (t = 0 ) = 1 } . 

Change of the integer variable —)■ rumax in this limit leads to 




sine {t) = lim 


mrnax-^-oo 777, 


E 

m=l 


777 — 1/2 
COS 1 - 1 


This signifies that if the integer r?7max is large enough, it retains all properties 
required to approximate the sine function that can be used for sampling (see 
sampling methodology in [23] for details). Consequently, we can generalize 
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the approximation ([^ of the complex error function for an arbitrary integer 
’^max as follows 


w{z) = W (z) ^ ^ 


+ {z + iz,/2) 


Im [ 2 ;] ^ 0. 


Cl-{z + iq/2f 
where the corresponding coefficients are rewritten as 

'vr (m — 1/2) [nh + q/2) 


(9) 


_ 1/2) A 

■^m r\ ^ 1 / J 


2m‘i.^h 




n=-N 


'^ma.xh 


Brr, = 


rrir, 


N 

E' 

n=-N 


PQg 1 7r(m-l/2) (n/i + ^/2) 




Cm. — 


TT (m — 1/2) 


2mmaxh 

Combining identity Q and approximation (|^ together at ?/ = 0 yields 
an exponential function approximation 


exp \ —x' 




+ (a^ + i'i/‘2‘) Bm ^ Am + {—X + k/2) B„ 


El L A - (a; + *‘^/2) C'i - {-X + iq/2y 


Figure 1 shows the difference £ (t) between the original exponential func¬ 
tion exp (—t^) and its approximation 


£ (t) = exp [—C) 


1 

2 


^max 

E 


m=l 


Am -f- (a; -f i<^/2) Bm ^ Am + {—x + i‘i/2y Bm 
- Cl^-{x + k/2f {-x + k/2f _ 


As we can see from this hgure, even with only rUmax = 16 summation 
terms the difference e (f) is very small and remains within the narrow range 
±5 X 10“^°. This conhrms a rapid convergence of the exponential function 
approximation that makes it suitable for numerical integration. Specihcally, 
this series approximation can be further used to replace the original exponen¬ 
tial function exp {—C) from the integrand in integral equation Q as follows 


K {x,y) « 


00 

/« -I "^max 

y 1 y- 

Am -t (t -t A/ 2 ) Bm 1 Am “t (~t + A/ 2 ) Bm 

27 r 7 y^ + {x- tf 

— 00 



( 10 ) 
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Fig. 1. The difference £ (t) between the original exponential function 
exp (— and its approximation at mmax = 16. 


Consider the series approximation (10) of the Voigt fnnction in more 
detail. The integrand in this integral is analytic everywhere over the entire 
complex plain except 2 + 4mmax isolated points 


{x iy, X + iy, Cm ^*5/2, Cm Cm T Cm T f?/2} ; 

TTl G {1, 2, 3, ... Tflmax} 


where singnlarities are observed. However, as we take a contonr integral only 
on the npper complex plane, for example as a semicircle Cccw with inhnite 
radius in counterclockwise (CCW) direction, the quantity of isolated points 
is reduced twice and becomes equal to 1 + 2mmax- 

Lastly, substituting the corresponding isolated points inside the domain 
enclosed by contour Cccw'- 


tr {x T iy, Cm T Cm T 5 rn G {1, 2, 3, ... nimax} 

into the Residue Theorem’s formula that for our specihc case is expressed in 
form 


^ f* 1H“ 2?72niax 

— (j) f{t)dt= Res [/ (t), tr], 

r=l 

'^ccw 


where / {t) is the integrand of integral (10), we hnd a new series approxima- 
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tion of the Voigt function 


K {x,y) 



'Cl-x^ + {y + ,l2f 

+ iBm {y + S‘/2) 

'Cl+x^ + {y + ,l2f 

[Cm + X — i {y + S‘/2)] [Cm “ 

-x + i{y + (,/2)] 

Cl,-{x + i{y + <,/2)f 


( 11 ) 


Since an algorithm involving complex numbers takes extra time, it would 
be very desirable to exclude them in computation. Thus, after some trivial 
rearrangements of the equation above, it can be represented in a simplihed 
form as the series approximation consisting of the real variables and constants 
only 

( \ A (/^m + y^ - X^) + '^raV + y^) 

^ ^iV) — 2_^ /Q2 _L OR i^.2 _ _L (rr.2 I 


+ 2 / 3 ™ (|/2 - a:2) + (3^2 ^ ^2)2 

K {x,y) K{x,y + q/ 2 ) , 


where 


2I™ 


a/tt (m - 1/2) 
2mL^^h 


N 

E 

n=-N 


/4 _ n'^h^ 

sm 


TT (m — 1/2) (nh + <^/2) 


'f^maxh 


o _r^-z _ / 7r(m- l/2) y 


and 


qm — f-Bm — 




N 

E 

n=-N 


cos 


71 {m — 1/2) (nh + ^/2) 

‘^maxh 


As the constants a™, /3™ and 7™ are independent of the input parameters x 
and y, the obtained series (12) is a rational approximation. 


3 Implementation 

Since the Voigt function is an even with respect to the parameter x and odd 
with respect to the parameter y. 


K (a;, - bl) = K {-x, - |r/|) = -K {x, |?/|), 
















it is sufficient to consider the values x and y only from the P* and 
quadrants in order to cover the entire complex plane. Consequently, in algo¬ 
rithmic implementation it is reasonable to take the second input parameter 
by modulus as |?/| and compute the Voigt function according to the scheme 


K{x,y>0)Ri k{x, \y\ + <^/2) 

K{x,y <0) ^ -K (x, \y\ + <^/2) 

Thus, if the parameter y is negative, we hrst take it by absolute value and, 
after computation, simply change the sign of the computed result to opposite. 
It should also be noted that taking the argument |?/| is advantageous in imple¬ 
mentation as it prevents computational overflow and enables an algorithmic 
stability (see Appendix A for details). 


The series approximation (12) alone covers the domain 0 < x < 40,000 
and 10’ <y < 10^, required in applications using the HITRAN molecular 
spectroscopic database [21]. In general, it provides accurate results while 
y ^ 10“®. However, this approximation may be used only to cover a smaller 
domain 0 ^ x ^ 15 and 10“® ^ ^ 15 that is considered most difficult for 

rapid and accurate computation. 

In our recent publication we have shown that the following approximation 
(see equation (|^ in ESI) can be effective for computation in the narrow 
domain 0 ^ x ^ 15 and 0 ^ y < 10“® along x-axis: 


K (x, y « 1) = Re [w (x, y « 1)] 


Re 


1 + 


le 


2F (x) — 


_ ^2ixy 


X 


or 


K {x,y « 1) 




cos(2xy) 



[y sine {2xy) 


F (x) sin {2xy)] , 


where 

X 


F(x) 



/ 


dt 


0 

is the Dawson’s integral. As argument x in the Dawson’s integral is real, its 
implementation is not problematic and several efficient approximations can 
be found in literature [2611271EH]. 
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When the input parameters x and y are large enough (say when the 
condition \x + iy\ > 15 is satished), many rational approximations become 
effective for accurate and rapid computation. For example, the Gauss-Hermit 
quadrature or the Taylor expansion can be effectively implemented (see for 
example [1] for details). 

A Matlab source code for computation of the Voigt/complex error func¬ 
tion that covers the entire complex plane can be accessed through Matlab 
Central, hie ID: [29]. This code has been developed by our research 

group and can be used for verihcation of the computed results. The domain 
divisions for computation of the Voigt function with complete coverage of 
the complex plane can be developed similarly. 

In order to demonstrate the computational efficiency of the series approx¬ 
imation (12), the comparison with the Weideman’s rational approximation 
has been made (see equation (38-1) and corresponding Matlab code in [T9]i. 
Such a choice is justihed since the Weideman’s approximation is one of the 
most rapid for computation of the Voigt/complex error function. The com¬ 
putational testing we performed by using a typical desktop computer shows 
that with same number of the summation terms rumax = 16 (default integer 
in Matlab code in [T9| is also 16), the algorithm based on series approxi¬ 


mation (12) is faster in computation than that of based in the Weideman’s 
rational approximation by factors about 2.2 and 2.7 for input arrays x and y 
consisting of 5 and 50 million elements, respectively (see the Matlab source 
code with implementation of the series approximation (12) in Appendix B). 
This is mainly because the Weideman’s rational approximation computes 
simultaneously both the real K [x, y) and imaginary L [x, y) parts, while 
the rational approximation (12) computes only the real part K {x, y) of the 
complex error function w{x,y). It should be noted that in most practical 
applications the imaginary part L {x,y) ([^ of the complex error function is 
not needed and simply ignored. Moreover, due to rapid convergence of the 
series approximation (12) we may decrease the number of the summation 
terms. In particular, at rumax = 12 the computational acceleration of the 
Voigt function can be further gained by about 30%. Therefore, the appli¬ 
cation of the series approximation (12) may be advantageous especially for 
intense computations with extended input arrays. 
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4 Error analysis 


In order to quantify accuracy of the series approximation (12), it is convenient 
to define the relative error as 


A = 


K (x, y) - Kref. (x, y) 


K-ref. {x,y) 


where Kref. (x, y) is the reference. The highly accurate reference values can 
be generated, for example, by using the Algorithm 680 [T8l l30] or recently 
published Algorithm 916 [HT] . 

]A of the relative error of the 
= 16. The domain required for coverage of 


Figures 2a and 2b show the logarithm log^gi 


series approximation (12) at m 


the HITRAN molecular spectroscopic database is 0 < a; < 40, 000 and 10’ -"< 
y < 10^ [3 |32] while the domain 0 ^ a; ^ 15 and 10“® ^ ^ 15 is the most 

difficult for accurate and rapid computation of the Voigt function. Therefore, 
we will consider the accuracy behavior within the HITRAN subdomain and 
narrow band domain 0 ^ a: ^ 15 fl 10“^ ^ ^ 15 and 0 ^ a: ^ 15 H 10“® ^ 

y ^ 10“^ separately as shown in Figs. 2a and 2b, respectively. 

As we can see from Fig. 2a, within the HITRAN subdomain the accu¬ 
racy of the series approximation is quite uniform and better than 10“^^ over 
most of this area. Although the accuracy deteriorates with decreasing y, it, 
nevertheless, remains high and better than 10“®. Another advantage is that 
the area where the accuracy deteriorates is relatively small. Particularly, the 
area where accuracy is worse than 10“^^ (yellow and red colors) is smaller 
than 2% of the domain’s total area. 

With randomly taken input parameters x and y, it is determined that 
the average accuracy over the domain of practical interest 0 < x < 40, 000 
and 10“^ < y < 10^ is 10“^^. Although the series approximation (12) can 


cover this domain accurately, it may be implemented only within domain 
0 ^ X ^ 15 and 10“® ^ ^ 15 that is the most difficult for accurate and 

rapid computation. 

In the narrow band shown in the Fig. 2b, the accuracy deteriorates 
further with decreasing y. However, it still remains high and better than 
10“®. In particular, the best and worst accuracies in the narrow band domain 
0 ^ X ^ 15 n 10“® ^ 1/ ^ 10“"^ exceed 10“^° (yellow color) and 10“® (dark 
red color), respectively. 

In modern applications requiring the HITRAN molecular spectroscopic 
data-base, the accuracy of the Voigt function should be 10“®. Therefore, we 
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Fig. 2. Logarithms of the relative error logj^o'^ 3) the 
HITRAN subdomain 0 ^ a; ^ 15 fl 10“"^ ^ y ^ 15 and b) for the 
narrow band domain 0 ^ a; ^ 15 fl 10“® ^ 1/ ^ lO”"^. The constants 
applied in computation are = 2.75, N = 23, rrtmax = 16 and 
h = 0.25. 


may reduce the integer rumax in the series approximation (12) from 16 to 12 
in order to gain computational acceleration. The number of the summation 
terms, determined by the integer rumax, is quite sensitive to the small param¬ 
eter value h. We have found empirically that at rumax = 12 the best accuracy 
can be achieved by taking h = 0.293. 

Figure 3a depicts the logarithm log^gA of the relative error of the series 
approximation (12) at rrimax = 12 in the HITRAN subdomain 0 ^ a; ^ 15 
and 10“^ ^ y ^ 15. One can see that in the HITRAN subdomain the 
accuracy is better than 10“®. 

Figure 3b illustrates the logarithm log^gA of the relative error of the series 
approximation (12) at rrimax = 12 in the narrow band domain 0 ^ a; ^ 15 
and 10“® ^ ^ 10“"^. Despite only 12 summation terms involved in the 

series approximation (12), the accuracy within this domain is better than 
10“®. For comparison, to achieve the same accuracy 10“® at y ^ 10“®, the 
Weideman’s approximation requires 32 summation terms (see Fig. 4 in [22] 
for details). Thus, we can see that the series approximation (12) may be 
useful and convenient in spectroscopic applications. 
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parameter x 


Fig. 3. Logarithms of the relative error 2 ) for the HITRAN 

subdomain 0 ^ a; ^ 15 fl 10“^ ^ y ^ 15 and b) for the narrow band 
domain 0 ^ a; ^ 15 fl 10“® ^ ^ 10“"^. The constants applied in 

computation are c = 2.75, N = 23, rrimax = 12 and h = 0.293. 


5 Conclusion 


A rational approximation for rapid and accurate computation of the Voigt 
function is presented. With only 16 summation terms, the proposed rational 
approximation provides average accuracy 10“^'^ in the domain of practical in¬ 
terest 0 < X < 40, 000 and 10“^ < y < 10^ that is needed for applications us¬ 
ing the HITRAN molecular spectroscopic database. The computational test 


shows that the algorithm based on series approximation (12) is more rapid in 


computation than that of based on the Weideman’s rational approximation 
by factor greater than 2. Algorithmic stability is achieved since the proposed 


series approximation (12) contains no poles at y ^ 0 and —00 < x < 00 . 
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Appendix A 


According to definition of the ^-function (12) we can write the following 
identity 


K{x,y + <,/2) = 


■^m 

Cl-x^ + {y + C2f 

+ iBm {y -t ^/2) 

Cm + x"^ + {y + <^/2) 

[Cm + X — i (2/ + <r/2)] [Cm “ 

-x + i{y + <;/2)] 

C^-{x + i{y + ^/2)f 


(A.l) 


and since 


[Cm + x-i{y + <;/2)] [Cm -x + i{y + <;/2)] \c^ - {x + i{y + f/2))^j = 


2 A-2) 

Pm + 2/3m (y{y + ?/ 2 )^ — + {y + ?/ 2 )^^ , 


where the right side of the identity (A. 2 ) is the denominator of k{x, 
y + q/2), it is sufficient to show that the poles do not exist on the right side 
of the identity (A.l) when both variables x and y are real such that ^ 0 in 
order to prove that K{x,y + <^/ 2 ) has no poles under same conditions. 


The proof is not difficult. Let us equate the left side of identity (A.2) to 


zero 


\Cm -\- X — i{y c/2)] [Cm — X i {y ‘^/2)] C‘^ — (x + z (y + ‘^/2)) — 0. (-^-3) 

and then solve this equation with respect to the variables x and y. Sup¬ 


pose now that the solutions in the equation (A.3) for real valued argu¬ 
ments X G (—oo, cx)) and y G [0, oo) exist. Solving the equation (A.3) 


with respect to x results in four possible solutions xi = —i{y + q/ 2 ) — Cm, 
X 2 = i{y + ?/ 2 ) Cm, X 3 = -i{y + ?/ 2 ) -F Cm and X 4 = i {y + ?/ 2 ) - Cm- 
Since the constants Cm, ^ are real valued and since ^ > 0, ^ 0, these 

solutions {xi, a; 2 , Xs,X 4 } must be always complex. However, the complex so¬ 
lutions {xi,X 2 ,X 3 ,Xi} contradict our initial assumption that x is real. Simi¬ 


larly, four possible solutions of equation (A.3) with respect to the variable y 
are yi = -i{x- Cm) - <^/ 2 , y 2 = i{x- Cm) - ?/ 2 , ys = -i{x + Cm) - <t /2 
and y 4 = i {x + Cm) — ‘^/2. Since the constants Cm, are real and positive, 
the solutions {|/i,I/ 2 , 1 / 3 , 2 / 4 } must be either complex at x 7 ^ Cm or negative 
and equal to —‘^/2 at x = Cm- However, the complex or negative solutions 
{yi,y 2 ,y 3 ,yi} contradict our initial assumption that y ^ 0. Due to these 
contradictions we must conclude that there are no poles in identity (A.3) 
under the conditions {x, ?/} G M such that y ^ 0- 
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The absence of the poles signihes that while the arguments is taken by 
absolute value as \y\, the function k (x, \y\ +?/2) will never encounter divi¬ 
sion to zero that leads to computational overflow. That is why taking the 
input parameter by absolute value as \y\ is advantageous since this approach 
provides stable performance of the algorithm. 


Appendix B 

function VF = voigtf(x,y,opt) 

y. This function file is a subroutine for computation of the Voigt function, 
y The input parameter y is used by absolute value according to the 
y procedure described in the article. The parameter opt is either 1 for 

y more accurate or 2 for more rapid computation. At y < 0 change the sign 

y to negative externally, out of the body of this function file. 

y 

y NOTE: This programi completely covers the domain 0 < x < 40,000 and 
y 10"-4 < y < 10~2 required for applications using the HITRAN molecular 
y spectroscopic database. However, it may be implemented only to cover the 
y smaller domain 0 <= x <= 15 and 10~-6 <= y <= 15 that is the most 

y difficult for rapid and accurate computation. See the article that 

y briefly describes how other domains can be covered. 

y 

y The code is written by Sanjar M. Abrarov and BrendcUi M. Quine, York 
y University, Canada, March 2015. 

if nargin == 2 
opt = 1; 

end 


if opt ~= 1 kk opt ~=2 

disp([’opt = ’,num2str(opt),’ cannot be assigned. Use either 1 or 2.’]) 
return 

end 

y Define array of coefficients as coeff = [alpha;beta;gamma]’ 
if opt == 1 
coeff = [ 

1.608290174437121e-001 3.855314219175531e-002 1.366578214428949e+000 

6.885967427017463e-001 3.469782797257978e-001 -5.742919588559361e-002 


15 



2.651151642675390e-001 
•2.050008245317253e-001 
•1.274551644219086e-001 
•1.134971805306579e-002 
4.201921570328543e-003 
8.084740485193432e-004 
1.946391440605860e-005 
•4.132639863292073e-006 
•2.656262492217795e-007 
•1.524188131553777e-009 
2.239681784892829e-010 
4.939143128687883e-012 
4.692078138494072e-015 
•2.512454984032184e-016 
]; 


9.638285547938826e-001 
1.889103967396010e+000 
3.122804517532180e+000 
4.664930205202391e+000 
6.515481030406647e+000 
8.674456993144942e+000 
1.114185809341728e+001 
1.391768433122366e+001 
1.700193570656409e+001 
2.039461221943855e+001 
2.409571386984707e+001 
2.810524065778962e+001 
3.242319258326621e+001 
3.704956964627684e+001 


-5.709602545656873e-001 
-2.011075414803758e-001 
1.069871368716704e-002 
1.468639542320982e-002 
1.816268776500938e-003 
-6.875907999947567e-005 
-2.327910355924500e-005 
-1.004011418729134e-006 
2.304990232059197e-008 
2.275276345355270e-009 
3.383885053101652e-011 
-4.398940326332977e-013 
-1.405511706545786e-014 
-3.954682293307548e-016 


mMax = 16; 16 summation terms 


elseif opt == 2 


coeff = [ 

2.307372754308023e-001 
7.760531995854886e-001 
4.235506885098250e-002 
-2.340509255269456e-001 
-4.557204758971222e-002 
5.043797125559205e-003 
1.180179737805654e-003 
1.754770213650354e-005 
-3.325020499631893e-006 
-9.375402319079375e-008 
8.034651067438904e-010 
3.355455275373310e-011 
]; 


4.989787261063716e-002 
4.490808534957343e-001 
1.247446815265929e+000 
2.444995757921221e+000 
4.041727681461610e+000 
6.037642585887094e+000 
8.432740471197681e+000 
1.122702133739336e+001 
1.442048518447414e+001 
1.801313201244001e+001 
2.200496182129099e+001 
2.639597461102705e+001 


1.464495070025765e+000 
-3.230894193031240e-001 
-5.397724160374686e-001 
-6.547649406082363e-002 
2.411056013969393e-002 
4.001198804719684e-003 
-5.387428751666454e-005 
-2.451992671326258e-005 
-5.400164289522879e-007 
1.771556420016014e-008 
4.940360170163906e-010 
5.674096644030151e-014 


mMax = 12; 12 summation terms 

end 


varsigma = 2.75; °L define the shift constant 
y = abs(y) + varsigma/2; 

arrl = y."2 - x.~2; define 1st repeating array 
arr2 = x."2 + y.~2; define 2nd repeating array 
arr3 = arr2.~2; Z define 3rd repeating array 

VF = 0; y. initiate VF 
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for m = l:mMax 

VF = VF + (coeff(m,1)*(coeff(m,2) + arrl) + ... 

coeff(m,3)*y.*(coeff(m,2) + arr2))./(coeff(m,2)~2 + ... 
2*coeff(m,2)*arrl + arr3); 

end 

end 
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